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ABSTRACT 

The  modelling  of  solid  propellant  ignition  is  investigated  with  the  aim  of  im¬ 
plementation  into  the  numerical  code  Casbar.  The  current  state  of  the  art  in 
solid  propellant  ignition  and  combustion  modelling  is  reviewed,  with  a  sim¬ 
plified  condensed  phase  ignition  model  chosen  as  the  suitable  candidate.  A 
number  of  methods  of  solving  the  mathematical  models  are  analysed,  with  re¬ 
sults  compared.  Based  on  these  results,  the  solution  by  integral  methods  to 
the  condensed  phase  model  was  chosen  for  implementation  into  Casbar. 
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Review  of  Solid  Propellant  Ignition  Models  Relative  to  the 
Interior  Ballistic  Modelling  of  Gun  Systems 

Executive  Summary 


Solid  propellant  ignition  (and  the  evolution  to  self-sustained  combustion)  is  a  highly 
complex  physicochemical  process,  involving  the  transition  of  a  stable  solid  propellant  state 
through  to  a  luminous  burning  resulting  from  the  application  of  heat  energy  [1],  For  solid 
rocket  and  gun  applications,  ignition  stimulus  is  usually  in  the  form  of  a  pyrotechnic  igniter 
emitting  heat  energy  and  hot  particles  into  the  propellant  bed.  The  action  of  simultaneous 
energy  sources  from  convection  of  the  hot  gases,  conduction  from  impingement  of  hot 
particles,  radiation  from  both  the  hot  igniter  gases  and  particles,  and  even  heat  from  atom 
recombination  and  vapour  condensation,  all  act  to  ignite  the  propellant  [2]. 

Interior  ballistic  modelling  is  used  in  a  wide  range  of  defence  applications,  and  forms 
a  key  analytical  tool  for  the  assessment  of  gun  and  rocket  propulsion  systems.  A  range 
of  phenomena  occurring  during  the  interior  ballistics  cycle  are  related  to  solid  propellant 
ignition  processes,  therefore  the  accurate  reproduction  of  ignition  phenomena  is  important. 
Australia’s  Defence  Science  and  Technology  Organisation  has  capability  in  performing 
gun  interior  ballistics  modelling  through  its  numerical  code,  Casbar.  Casbar  solves  the 
governing  equations  for  the  transient  flow  of  chemically  reacting  gas  and  particulates  within 
a  finite  volume  discretisation  of  the  computational  domain  [3]. 

Currently  in  Casbar  propellant  ignition  is  modelled  using  a  simple  go/no-go  condition, 
whereby  if  the  gas  surrounding  the  propellant  is  above  the  defined  propellant  ignition 
temperature,  the  grain  will  combust.  This  ignition  criterion  does  not  account  for  finite- 
rate  grain  heating  and  the  experimentally  observed  ignition  delay  of  solid  propellant  grains, 
and  thus  generally  predicts  an  unrealistically  fast  propellant  ignition. 

This  report  describes  the  solid  propellant  ignition  and  combustion  phenomena,  and 
investigates  a  number  of  ignition  models  suitable  for  implementation  into  the  Casbar  code. 
The  three  main  areas  of  ignition  models  (solid-phase,  heterogeneous  and  gas-phase  reac¬ 
tions  models)  encompass  a  broad  range  of  complexity  and  numerical  efficiency.  It  is  desir¬ 
able  to  choose  a  model  that  is  accurate,  while  being  flexible.  Ultimately,  the  solid-phase 
ignition  model  was  chosen  for  implementation  in  Casbar. 

A  number  of  numerical  techniques  for  solution  of  the  solid-phase  ignition  model  were 
reviewed.  In  numerical  modelling  it  is  required  that  the  solution  of  the  model  be  adequately 
robust,  while  reducing  the  impact  on  the  simulation  time.  In  all  of  the  investigated  heat 
flux  scenarios,  the  integral  method  was  able  to  adequately  approximate  the  final  ignition 
time  to  within  an  acceptable  level  of  accuracy  (in  comparison  to  the  finite  difference  ap¬ 
proximation).  Situations  involving  highly  variable  heat  fluxes  were  employed  to  test  the 
applicability  of  the  integral  method.  These  scenarios  were  constructed  to  accentuate  the 
variability  of  the  heat  flux,  and  situations  like  this  are  not  expected  in  typical  interior  bal¬ 
listic  simulations.  The  integral  method  is  therefore  considered  an  appropriate  candidate 
for  implementation  in  Casbar. 
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1  Introduction 

The  accurate  reproduction  of  ignition  phenomena  is  an  important  aspect  of  interior  ballistic 
modelling.  A  range  of  phenomena  occurring  during  the  interior  ballistics  cycle  are  related 
to  solid  propellant  ignition  processes.  Australia’s  Defence  Science  and  Technology  Organ¬ 
isation  currently  has  capability  in  performing  gun  interior  ballistics  modelling  through  its 
code,  Casbar.  Casbar  solves  the  governing  equations  for  the  transient  flow  of  chemically 
reacting  gas  and  particulates  within  a  finite  volume  discretisation  of  the  computational 
domain  [3].  Currently  propellant  ignition  is  modelled  using  a  simple  go/no-go  condition, 
whereby  if  the  gas  surrounding  the  propellant  is  above  the  defined  propellant  ignition  tem¬ 
perature,  the  grain  will  combust.  This  ignition  criterion  does  not  account  for  finite-rate 
grain  heating  and  the  experimentally  observed  ignition  delay  of  solid  propellant  grains,  and 
thus  generally  predicts  an  unrealistically  fast  propellant  ignition.  This  report  describes  the 
solid  propellant  ignition  and  combustion  phenomena,  and  investigates  a  number  of  ignition 
models  suitable  for  implementation  into  the  Casbar  code. 


1.1  Solid  propellant  ignition  and  combustion 

Solid  propellant  ignition  (and  the  evolution  to  self-sustained  combustion)  is  a  highly  com¬ 
plex  physicochemical  process,  involving  the  transition  of  a  stable  solid  propellant  state 
through  to  a  luminous  burning  resulting  from  the  application  of  heat  energy  [1].  For  solid 
rocket  and  gun  applications,  ignition  stimulus  is  usually  in  the  form  of  a  pyrotechnic  igniter 
emitting  heat  energy  and  hot  particles  into  the  propellant  bed.  The  action  of  simultane¬ 
ous  energy  sources  from  convection  of  the  hot  gases,  conduction  from  impingement  of  hot 
particles,  radiation  from  both  the  hot  igniter  gases  and  particles,  and  even  heat  from  atom 
recombination  and  vapour  condensation,  all  act  to  ignite  the  propellant  [2]. 

Gun  and  rocket  propellants  traditionally  comprise  several  chemical  ingredients  com¬ 
bined  either  heterogeneously  or  homogeneously  to  form  the  solid  [4].  The  chemical  ingredi¬ 
ents,  which  can  consist  of  oxidizer,  fuel,  binder,  plasticiser,  curing  agent,  stabilizer,  bonding 
agent,  burning  rate  catalyst,  anti-aging  catalyst,  opacifier,  flame  suppressant,  combustion 
instability  suppressant  and  cross-linking  agent  [5]  are  either  connected  on  a  microscopic 
level  (homogeneous)  or  mixed  on  a  macroscopic  level  (heterogeneous)  to  form  the  desired 
physical  and  chemical  properties  of  the  propellant.  The  ignition  and  combustion  processes 
for  homogeneous  and  heterogeneous  propellants  possess  many  subtle  differences  between 
formulations  and  ingredients,  however  the  underlying  physical  mechanisms  are  similar. 
Traditionally,  homogeneous  propellants  (such  as  single  base  nitrocellulose  or  double  base 
nitrocellulose  and  nitroglycerine)  are  employed  in  gun  propellants,  while  heterogeneous 
(such  as  AP,  HMX  or  RDX)  propellants  are  more  often  employed  in  solid  rocket  propel¬ 
lants.  Composite  propellants  are  finding  more  use  in  gun  propellant  applications  in  recent 
times  [6],  however  AP  based  propellants  are  not  used  in  gun  applications  due  to  the  large 
amount  of  hydrogen  chloride  (HC1)  produced  in  the  combustion  products,  which  acts  to 
accelerate  gun  barrel  erosion  [7]. 

Many  complex  processes  occur  between  the  initial  application  of  heat  energy  and  steady 
state  combustion;  solid  propellants  are  converted  to  gas  products  through  oxidation  at  ex¬ 
tremely  high  pressures  and  temperatures  in  the  order  of  a  few  thousand  Kelvin  can  be 
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reached.  All  of  this  occurs  over  a  fraction  of  a  second,  resulting  in  an  extremely  hostile 
environment  that  is  still  not  completely  understood,  despite  decades  of  research  [8].  Sig¬ 
nificant  condensed  phase  reactions,  coupled  with  a  large  number  of  gas  phase  reactions, 
compound  the  difficulties  of  forming  a  suitable  model  of  solid  propellant  ignition  and  com¬ 
bustion. 

When  a  propellant  is  subjected  to  sufficient  heat  energy  such  that  self-sustained  com¬ 
bustion  is  obtained,  three  significant  reaction  regions  develop  defined  by  the  thermal  state 
of  the  propellants  (or  its  constituents),  as  shown  in  Figure  1.  The  solid  (or  condensed) 
phase  region  is  the  furthest  from  the  source  of  heat  energy,  and  consists  of  the  propellant  in 
its  initial  solid  state.  Upon  the  external  application  of  heat  the  surface  temperature  begins 
to  increase  in  conjunction  with  heat  being  conducted  away  from  the  surface  into  the  solid 
core.  The  surface  temperature  continues  to  increase  until  the  point  of  phase  change  is 
reached.  Solid-phase  exothermic  reactions  may  occur  in  some  propellant  ingredients  (such 
as  AP  [9]  or  ADN  [4])  leading  to  thermal  degradation  within  the  solid  phase,  however  these 
are  typically  insignificant  and  often  neglected.  The  second  region  can  generally  be  consid¬ 
ered  to  consist  of  either  a  two-phase  mixture  of  melted  liquid  propellant  and  evaporated 
propellant  gas  or  an  infinitesimally  small  sublimation  interface.  Some  propellants  exhibit 
melting  behaviour  while  others  do  not  [10],  in  which  the  melt  layer  can  be  formulated  with 
sublimation  or  pyrolysis  mechanisms  for  transition  to  the  gas  phase  [11].  Although  the 
melt  line  is  often  represented  as  a  sharp  discontinuity  in  theoretical  models,  in  reality  it  is 
a  blurred  region  consisting  of  a  slurry  of  solid  and  liquid  propellant.  There  is  no  definite 
distinction  between  regions,  and  complex  reactions  such  as  vaporization,  condensation,  de¬ 
composition  and  oxidation  can  all  occur.  This  layer  is  often  referred  to  as  the  foam  layer 
due  to  its  ‘frothy’  nature  [7].  The  final  region  in  the  combustion  process  is  the  gas-phase 
layer,  a  region  where  the  evaporated  species  react  and  decompose  into  other  species. 

The  gas  layer  itself  has  transient,  pressure  dependent  behaviour  that  can  consist  of 
multiple  sub-layers.  At  low  burning  pressures  (below  approximately  1  MPa)  there  is  no 
visible  flame  above  the  burning  surface  of  the  propellant.  As  the  pressure  increases,  a 
■weak  flame  is  sometimes  visible  above  the  burning  surface.  As  the  pressure  increases 
further,  the  flame  moves  closer  to  the  propellant  surface,  until  pressures  of  around  lOMPa, 
where  it  appears  attached  to  the  propellant  burning  surface  [10].  The  region  between  the 
burning  surface  and  visible  flame  can  consist  of  two  regions,  labelled  the  fizz  zone  and 
the  dark  zone.  Within  both  these  regions,  the  gas  temperature  is  relatively  constant  due 
to  a  slow  reduction  of  NO  to  N2  [10].  Rapid  reactions  occur  in  the  fizz  zone  just  above 
the  propellant  surface.  Within  the  dark  zone  the  oxidised  products  released  from  the  fizz 
zone  react  relatively  slowly,  and  the  zone  is  only  present  if  the  pressure  and  temperature 
are  sufficient  [7].  Within  the  flame  zone,  the  oxidation  reactions  release  large  amounts 
of  energy  with  the  final  temperature  approaching  the  adiabatic  flame  temperature  of  the 
propellant.  The  heat  feedback  from  the  flame  in  the  form  of  radiation,  convection  and 
conduction  sustains  the  reactions  until  the  solid  propellant  is  consumed.  Flame  stand-off 
plays  a  critical  role  in  the  amount  of  heat  feedback  to  the  propellant. 
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Figure  1:  Solid  propellant  combustion  schematic  (adpated  from  [4j). 
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2  Current  state  of  solid  propellant  ignition  and 

combustion  modelling 

Propellant  ignition  modelling  has  been  an  area  of  research  since  the  1950s,  with  many 
different  theoretical  models  developed  [5].  Propellant  ignition  models  can  be  c.atagorized 
into  three  main  types:  solid-phase  (reactive  solid),  heterogeneous,  and  gas  phase  reaction 
models  [4]. 

Solid  phase  reaction  models  were  the  first  developed,  and  are  generally  the  least  com¬ 
plicated  due  to  the  extensive  simplifications  employed  [12].  The  internal  energy  of  the 
propellant,  often  coupled  to  a  single-step  Arrhenius  reaction  within  the  condensed  phase, 
is  considered  the  dominant  driving  force  of  ignition.  The  effects  of  gas  phase  reactions  are 
assumed  to  not  influence  the  time  to  ignition,  and  therefore  neglected. 

Heterogeneous  ignition  models  were  the  next  form  of  ignition  model  developed,  due  to 
experimental  observations  of  ignition  from  introduced  hypergolic  oxidizer  gases  [13].  These 
models  allow  for  observed  behaviour  of  composite  propellants  by  including  the  diffusion 
of  fuel  and  oxidizer  species  at  the  propellant  surface.  The  complexity  of  the  S3^stem  is  in¬ 
creased  over  the  simple  condensed  phase  models  through  the  tracking  of  species  concentra¬ 
tion  at  the  solid-gas  interface.  Heterogeneous  ignition  models  are  based  on  the  assumption 
that  surface  interface  reactions  are  the  driving  force  for  ignition.  The  gas  phase  is  assumed 
to  consist  of  products  from  the  combustion  process. 

The  third  type  of  solid  propellant  ignition/combustion  model  is  the  gas  phase  model. 
Gas  phase  models  increase  the  complexity  over  heterogeneous  models  by  incorporating 
the  reactions  and  species  concentration  within  the  gas  phase.  The  gas  phase  reactions  are 
considered  to  be  the  driving  factor  towards  ignition,  with  the  condensed  phase  contributions 
to  species  decomposition  generally  neglected.  The  gas  phase  models  have  the  ability  to 
treat  the  entire  ignition  transient  from  heat  addition  to  the  stable  solid  right  through  to 
pressure-dependent  steady-state  combustion  of  the  propellant.  Some  models  may  be  used 
to  predict  the  burn  rate  characteristics  of  composite  propellants,  however  these  are  still 
in  their  infancy  and  rely  heavily  on  empirical  results  [8].  This  section  gives  an  overview 
of  the  research  into  solid  propellant  ignition  and  combustion  models  from  the  early  solid 
phase  models,  through  to  the  current  multi-species  gas  phase  models. 


2.1  Solid-phase  (reactive  solid)  ignition  models 

Solid  phase  ignition  models  began  as  an  extension  of  the  solid  phase  thermal  explosion 
theory.  Hicks  [12]  developed  a  physical  model  of  the  homogeneous  solid  propellant  ignition 
process,  with  no  diffusion  or  consumption  of  propellant,  and  no  chain  reactions.  The  model 
is  assumed  to  be  one-dimensional,  with  the  heat  penetration  only  occurring  over  a  fraction 
of  the  solid  propellant  depth  for  rapid  application  of  heat.  This  assumption  has  been 
validated,  with  the  thermal  penetration  depth  observed  to  be  in  the  order  of  100 /xm  at 
solid  rocket  motor  operating  pressures  [4].  Although  solid  rocket  motor  operating  pressures 
are  significantly  lower  than  t3rpical  gun  systems,  during  the  ignition  phases  the  pressure 
within  the  gun  chamber  is  often  well  below  solid  rocket  motor  operating  pressures.  The 
one-dimensional  heat  conduction  equation  states  that  the  rate  of  energy  accumulation 
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within  the  solid  is  equal  to  the  rate  of  energy  conduction  plus  the  energy  production  due 
to  solid  phase  exothermic  reactions.  This  is  therefore  expressed  in  partial  differential  form 
as 


dT  ,  d2T  ■ 

=kw+Q 


(i) 


Q  is  the  rate  of  heat  evolution  per  unit  volume  due  to  solid  phase  reactions,  assumed 
to  be  of  zero  order  and  independent  of  species  concentration,  giving 


Q  =  Qrfe~EA/RT  (2) 

where  Qr  is  the  heat  of  reaction  per  unit  volume  (^/m3).  Boundary  conditions  were 
formulated  based  on  the  assumption  that  all  important  physical  processes  occur  within 
the  solid  propellant.  It  is  assumed  that  the  thermal  wave  penetration  is  small,  such  that 
the  temperature  of  the  propellant  at  a  relatively  large  distance  into  the  grain  is  always 
equal  to  the  initial  temperature  of  the  propellant,  To.  This  assumption  holds  true  if  the 
time  between  heating  and  ignition  is  less  than  the  time  for  the  thermal  wave  to  penetrate 
through  the  solid.  This  is  valid  for  most  gun  applications,  where  high  heating  rates  are 
applied  over  short  time  scales.  If  this  is  not  the  case,  the  boundary  condition  of  zero 
temperature  gradient  at  the  propellant  centre  can  be  applied.  At  the  solid/gas  boundary, 
the  heat  flux  is  taken  as  the  heat  feedback  from  the  gas  to  the  solid  phase,  such  that 

dT 

£^U=o  =  q(t)  (3) 

The  resulting  system  of  equations  was  solved  by  neglecting  solid  phase  reactions  (Q  =  0) 
and  transforming  the  variables  to  non-dimensional  form. 

Price  et  ol.  [14]  considered  the  ignition  of  solid  propellant  under  the  influence  of 
a  laser,  extending  the  model  to  accommodate  the  depletion  of  the  solid  material  at  a 
surface  regression  rate  r.  Optical  absorption  of  the  laser  energy  due  to  transparency  of 
the  propellant  within  the  energy  equation  for  the  solid  phase  was  also  included.  The  solid 
phase  reaction  partial  differential  equation  (PDE)  was  therefore  extended  to 


U1  U  1  ,U1  n.  _a  ■ 

pc"u^kd^  +  pvai  +  l3qe  +Q  (4) 

where  (3  is  the  absorptivity  of  the  laser  energy  into  the  solid  propellant  (m-1)  and  q 
is  the  energy  flux  per  unit  area  from  the  laser  ignition  source  ( w/m2-s ).  The  boundary 
conditions  are  similar  to  the  previous  model,  however  an  updated  heat  balance  at  the 
moving  boundary  is  introduced, 


Qi  +  prQs  =  k—  (5) 

where  g*  is  the  rate  of  heat  transferred  to  the  interface  and  Qs  is  the  heat  source  as¬ 
sociated  with  the  interface  (which  includes  both  chemical  reactions  and  phase  change). 
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The  boundary  condition  at  x  =  oo  can  be  either  taken  as  a  fixed  temperature  (as  per  the 
original  model)  or,  if  condensed  phase  reactions  are  present,  as  zero  gradient.  Although  an 
increase  in  the  realism  of  the  model  over  the  original  condensed  phase  model,  considerable 
simplifications  were  required  to  form  a  tractable  solution  [14].  The  thermophysical  proper¬ 
ties  (such  as  p.  cv,  k,  (3  and  Qs )  are  assumed  to  be  constant  and  therefore  do  not  vary  with 
time,  temperature,  pressure  or  propellant  composition  uniformity  and  optical  properties. 
The  assumption  of  single  step,  global  Arrhenius  rate  equation  for  solid  phase  reactions  is 
also  a  substantial  simplification  of  the  highly  complex  physicochemical  process  occurring. 
Even  to  the  current  day,  the  solid  phase  reactions  are  still  not  well  understood  [4]. 

Baer  and  Ryan  [15]  have  further  extended  the  solid  phase  analysis  to  include  the  tran¬ 
sition  from  ignition  to  steady  state  regression,  as  opposed  to  a  simple  go/no-go  condition. 
A  relationship  for  the  transition  is  produced  bj^  the  consideration  of  an  energy  balance 
at  the  surface  of  the  solid  propellant.  In  their  ignition  model,  the  regression  of  the  solid 
propellant  surface  under  a  known  heat  flux  is  taken  into  consideration,  and  is  a  function 
of  surface  temperature.  A  one-dimensional  heat  conduction  equation  neglecting  condensed 
phase  reactions  is  produced  for  dimensionless  variables; 

<9T  ,dT  d2T 

— -  —  r— —  =  — —  (6) 

dt  dx  Qx  2 

The  heat  feedback  q  to  the  solid  surface  is  assumed  to  occur  from  a  single  chemical 
reaction  at  (or  near)  the  propellant  surface.  This  feedback  is  related  to  the  heat  flux  in 
the  solid  and  the  regression  rate,  given  by 


q  =  qs  +  r-Q  (7) 

where  qs  is  the  dimensionless  heat  flux  just  below  the  solid  surface  and  r  is  the  dimen¬ 
sionless  rate  of  gasification.  Under  steady-state  regression,  the  dimensionless  heat  flux  just 
below  the  surface  is  related  to  the  difference  between  the  steady-state  regression  temper¬ 
ature  and  the  solid  initial  temperature,  multiplied  by  the  linear  regression  rate.  For  large 
values  of  the  steady  state  regression  temperature,  the  heat  feed  back  to  the  solid  surface 
is  approximately  related  to  pressure  b3^  an  exponential  factor.  During  the  early  ignition 
stages  (when  the  surface  temperature  is  much  lower  than  during  steady-state  regression), 
the  surface  temperature  is  assumed  to  be  dependent  on  the  rate  of  generation  of  the  reac¬ 
tive  species  from  the  solid  phase.  To  approximate  the  transition  from  a  kineticallv  limited 
condition  to  a  diffusion-limited  case,  the  following  relationship  is  established 


1  _  1  1 

q  qapn  e-1/Ts 


(8) 


where  qa  is  the  dimensionless  heat  feedback  at  1  atm  and  n  is  the  pressure  exponent  for 
steady  state  feedback  flux.  The  heat  flux  is  then  used  to  calculate  the  solid  temperature 
profile  through  a  Neumann  boundary  condition.  This  model  is  not  justified  theoretically 
by  the  authors,  it  is  employed  because  it  simply  produces  the  correct  as3rmptotes  as  the 
surface  temperature  is  varied,  while  providing  a  smooth  transition  between  the  two  rates. 
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Perez  et  al.  [16]  established  a  model  for  the  initial  pressure  rise  to  steady  state  combus¬ 
tion  (starting  transient)  of  solid  propellant  rocket  motors  with  high  internal  gas  velocities, 
with  comparison  of  their  model  with  experimental  data.  The  one  dimensional  heat  conduc¬ 
tion  equation  was  solved  for  the  solid  propellant  assuming  the  thermal  wave  penetration 
depth  was  negligible.  The  heat  transfer  to  the  propellant  surface  was  calculated  based 
on  a  local  convective  heat  transfer  coefficient.  A  third-degree  polynomial  solution  for  the 
propellant  surface  temperature  was  obtained  through  an  integral  method  [17],  which  was 
then  employed  to  calculate  the  surface  temperature.  The  propellant  was  considered  ignited 
when  the  surface  temperature  was  above  a  specific,  empirically  derived  value. 


2.2  Heterogeneous  ignition  models 

Heterogeneous  reaction  models  "were  born  from  the  experimental  observations  of  the  igni¬ 
tion  of  solid  propellants  subject  to  highly  energetic  oxidising  gas  environments  via  surface 
reactions  at  the  solid/gas  interface  [18].  This  theory  is  centered  around  exothermic  surface 
reactions  supplying  the  energy  for  stable  combustion,  and  should  not  be  confused  wdth  only 
applying  to  the  heterogeneous  nature  of  composite  propellants  [19].  In  these  models,  the 
rate  driving  reaction  was  assumed  to  occur  at  the  interface  of  the  solid  propellant  and  the 
oxidising  gas.  Diffusion  between  the  decomposed  solid  propellant  gasses  and  the  oxidizer 
provided  the  fuel  to  initiate  and  sustain  combustion.  Williams  [20]  analyzes  the  ignition 
of  a  solid  propellant  through  heterogeneous  reactions.  Similar  to  the  homogeneous  solid 
propellant  models,  a  one-dimensional  time  dependent  system  is  produced  that  involves  the 
reaction  between  the  solid  propellant  and  an  oxidising  gas  at  the  solid/gas  interface.  A 
known  heat  flux  is  introduced,  which  is  absorbed  at  the  solid/gas  interface.  Heat  con¬ 
duction  is  allowed  to  occur  in  both  the  solid  and  gas  phases,  however  no  phase  change  is 
accounted  for.  Condensed  phase  reactions  are  also  neglected,  however  a  single  step  Ar¬ 
rhenius  reaction  at  the  interface  is  present.  Ignition  is  said  to  occur  when  the  interface 
temperature  reaches  a  certain  value.  A  set  of  conservation  equations  is  produced  for  the 
solid,  gas  and  interface  regions.  For  the  solid  and  gas  phases,  the  conservation  of  energy 
equation  for  both  phases  is 


dh  _  d  (a  ■  dh/o^) 

dt  ~  dip  U 

where  a  =  p-  k/cp  and  ip  is  the  stream  function  ip  =  pdx.  The  conservation  of  energy 
condition  at  the  solid/gas  interface  (ip  =  0)  is  written  as 


dh* 


dha 


dhi 


as  dip  dip  Qi  dip  +q 


with  the  Arrhenius  rate  equation  at  the  interface  given  as 


(10) 


=  QBf  (Do  +  SQ)  e_TA/(T°+hg/cvg)  (11) 

where  B  is  a  pre-exponential  rate  factor,  Q  the  constant  heat  released  in  the  surface 
reaction  (per  unit  mass),  and  /  is  a  function  that  expresses  the  dependence  of  the  reaction 
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on  the  concentration  of  the  oxidizing  agent  in  the  gas  phase.  Williams  [20]  approximates 
this  as  f  (Y)  =  Yn ,  where  n  is  the  order  of  the  reaction.  A  Laplace  transform  is  then  used 
to  solve  the  PDEs. 


2.3  Gas-phase  reaction  ignition  models 

Gas  phase  models  generally  consider  all  (or  most)  of  the  physical  mechanisms  occurring 
during  propellant  combustion,  thereby  forming  a  complete  mathematical  description  of 
the  process.  Beckstead  et  al.  [4]  have  produced  a  very  detailed  review  of  the  gas  phase 
reaction  model,  with  a  summary  of  the  physical  processes  reproduced  here.  As  with  all  of 
the  models,  simplifications  may  be  included  to  reduce  the  complexity  of  the  system  to  be 
solved. 

The  solid  phase  region  is  modelled  in  a  similar  fashion  to  that  described  within  Section 
2.1,  with  the  conservation  of  energy  considered.  In  the  case  of  propellant  mixtures,  the 
properties  are  taken  as  mass  fraction  averages  of  the  constituents.  The  solid-phase  exother¬ 
mic  reactions  are  not  considered  to  contribute  significantly  to  the  ignition  transient,  and 
hence  are  neglected. 

External  to  the  solid  phase  region  is  the  two-phase  subsurface  region,  which  consists  of 
melted  liquid  propellant  and  evaporated  propellant  gas.  Within  this  region  there  are  a  large 
number  of  complicated  physicochemical  processes  occurring,  which  significantly  increase 
the  complexity  of  the  system.  Thermal  decomposition,  evaporation,  bubble  formation,  gas- 
phase  reactions  within  bubbles  and  transport  of  mass  and  energy  between  phases  are  all 
known  to  occur.  One  technique  of  simplifying  the  highly  complex  fluid  dynamic  interactions 
within  the  subsurface  region  is  to  employ  a  spatial  averaging  technique  [4].  A  fractional- 
voidage  0/  is  used  to  define  the  cross  sectional  areas  assumed  by  the  gas  phase  and  the 
condensed  phase,  assuming  the  bubbles  are  small  and  well  dispersed.  The  area  given  by 
the  gas  bubbles  is  therefore 


A9  =  (t)fA 

where  Ag  is  the  cross  sectional  area  of  the  propellant  gas  bubbles  and  A  is  the  cross 
sectional  area  of  the  propellant  sample.  The  integral  form  of  the  conservation  equations 
are  then  combined  using  the  fractional  voidage.  Species,  mass  and  energy  conservation 
equations  are  developed  for  the  liquid  and  gas  phases.  Within  the  subsurface  region  both 
evaporation  and  condensation  are  occurring  simultaneously.  The  net  difference  between 
the  condensation  and  evaporation  rates  is  determined  through  empirical  means  [21].  The 
mass  conversion  rate  due  to  evaporation  is  then  determined  by  multiplying  the  net  mass 
conversion  rate  by  the  specific  surface  area,  an  empirically  derived  relationship  between 
the  number  of  bubbles  and  the  fractional  voidage. 

The  gas  phase  analysis  is  centred  around  the  conservation  of  mass,  energy  and  species 
transport  for  a  multi-component  chemically  reacting  system  [4].  The  system  of  equations 
allows  for  the  accommodation  of  finite-rate  chemical  kinetics  and  variable  thermophysical 
properties.  A  multi-phase  treatment  similar  to  the  subsurface  region  can  be  employed  when 
the  gas  phase  contains  dispersed  condensed  phase  species.  All  thermophysical  properties 
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are  mass-averaged,  and  mass  diffusion  velocity  U)  of  an  individual  specie  arises  from  both 
concentration  and  temperature  gradients,  given  by 


Vi  =  - A 


1  dXj 
Xi  dx 


+  Dt 


DTz  1  dTg 
Xi  T  dx 


To  close  the  system  of  equations  the  ideal  gas  law  for  a  multi-component  system  is  used, 
coupled  "with  suitable  boundary  conditions  to  describe  the  energy  and  mass  conversion 
occurring  at  the  region  interfaces. 

As  the  gas  phase  reaction  model  is  the  most  detailed  of  the  models  examined,  with 
both  the  solid  and  gas  phase  considered,  a  simplified  formulation  is  highly  desirable.  The 
solid  phase  conservation  of  energy  equation  is  formulated  considering  surface  regression; 


dTs  .  dTs  ,  d2Ts 

Asc.s— +  Asrc.s— 

however  the  boundary  condition  at  the  propellant  surface  is  updated  to  include  the 
energy  associated  with  the  change  of  phase  from  solid  to  gas.  This  gives  the  boundary 
condition  at  x  =  0  as 


~dx  =«(t)+^-9 

where  is  the  heat  associated  with  phase  change.  The  sub-surface  region  is  assumed 
to  be  infinitely  thin  and  only  occur  in  the  interface  between  the  solid  and  gas  regions. 
Conservation  of  mass  within  the  solid  phase  (assuming  a  single  component  system)  is 
given  b3r 


ms  =  psf 

which  simply  states  that  the  rate  of  change  of  mass  within  the  solid  phase  is  due  to  the 
regression  of  the  propellant  surface,  multiplied  by  the  densit3'.  A  mass  balance  between 
the  solid  and  gas  phase  yields  the  conservation  of  mass  equation  for  the  gas  phase, 


Wig  =  PgUg  =  Tfls  =  psfS 

Hence,  the  rate  of  mass  addition  into  the  gas  phase  is  equal  to  the  mass  loss  from  the 
solid  phase  due  to  evaporation  or  sublimation.  At  the  propellant  surface,  an  equation  is 
required  that  describes  the  phase  change  process  to  determine  the  rate  of  regression,  r.  A 
number  of  models  have  been  produced  that  give  an  approximation  to  the  mass  transfer 
between  the  phases  covering  a  broad  range  of  complexity.  The  decomposition  of  the  solid 
has  been  approximated  through  the  inclusion  of  a  pressure  dependence  on  surface  regression 
[22],  The  rate  of  regression  of  the  solid  propellant  surface  is  given  by 


-ea/R 


Tsurf  ^  Tign 

Tsurf  Tign 
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where  Ea  is  the  activation  energy  for  the  pyrolysis  equation.  This  model  provides 
a  smooth  transition  from  the  non-combusting  decomposition  to  the  steady  state  burning 
condition. 

The  conservation  of  energy  equation  for  the  gas  phase  is  similar  to  the  solid  phase,  and 
is  given  by 


3Tg  .  dTg 

+  PcTC,-^ 


d2T 

l.  ° 

9  dx2 


+  qgwg  (x) 


w'here  qg  represents  the  heat  released  and  wg  is  the  mass  of  reactant  produced  by  the  gas 
chemical  reactions.  Single  or  multiple  specie  kinetic  mechanisms  can  be  considered.  Casbar 
currently  has  the  provision  for  multi-component  gas  phase  reactions,  allowing  complex  gas 
kinetics  of  propellant  combustion  to  be  modelled.  Arrhenius  rate  equations  are  used  for 
the  intermediate  reaction  steps.  The  boundary  condition  applied  to  the  gas  region  is 


Tx=o  —  Tsurf 

where  Tsurj  is  the  propellant  surface  temperature.  The  gas  temperature,  and  therefore 
resulting  heat  feedback  to  the  propellant  surface,  is  determined  through  the  chemical  kinetic 
model  and  the  conservation  equations  for  the  gas  phase.  In  this  model,  propellant  ignition 
is  assumed  to  have  occurred  when  the  heat  feedback  to  the  propellant  surface  is  able  to 
sustain  the  combustion  without  external  influence  however  this  definition  is  somewhat 
arbitrary  as  there  is  no  definitive  “ignition”  point. 


2.4  Defining  the  point  of  ignition 

The  establishment  of  deflagration  within  solid  propellant  is  difficult  to  define,  making  the 
specification  of  the  point  of  ignition  somewhat  vague.  Ignition  can  be  considered  to  have 
occurred  when,  upon  the  removal  of  the  ignition  stimuli,  the  propellant  is  able  to  achieve 
self-sustained  combustion  with  no  further  application  of  energy  [14].  From  a  modelling 
perspective,  it  is  only  important  to  have  a  distinct  ignition  criteria  if  the  ignition  and 
combustion  models  are  separate  (for  example  if  \rieille’s  lawT  is  employed  for  pressure- 
dependant  burning).  For  gas-phase  reaction  models  defining  an  ignition  point  is  somewhat 
arbitrary,  as  they  are  capable  of  modelling  the  transition  from  a  solid  propellant  at  ambient 
conditions  right  through  to  steady  state  combustion. 

Just  as  there  is  a  wide  range  of  theories  regarding  propellant  ignition,  there  are  also 
many  different  methods  for  determining  whether  propellant  ignition  has  occurred  in  both 
experimental  and  numerical  investigations  [19].  In  experimental  investigations,  often  pro¬ 
pellant  ignition  is  established  by  light  emission,  pressure  and/or  temperature  rise,  be¬ 
haviour  of  the  propellant  after  the  removal  of  ignition  stimuli  or  examination  of  the  pro¬ 
pellant  grain  after  quenching.  In  mathematical  investigations,  typically  ignition  is  defined 
to  have  occurred  when  a  certain  surface  temperature  has  been  reached  (the  ignition  tem¬ 
perature),  the  rate  of  temperature  rise  is  above  a  specified  value,  the  rate  of  heat  generation 
is  greater  than  the  rate  of  heat  loss  to  the  surroundings  or  b3^  observation  of  the  behaviour 
of  the  mathematical  model  after  the  artificial  ignition  stimuli  has  been  removed.  It  is 
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important  to  select  an  ignition  criterion  not  only  relevant  to  the  process  being  modelled, 
but  also  the  ignition  model  itself. 


2.5  Ignition  model  summary 

The  Casbar  code  was  developed  for  the  investigation  of  the  interior  ballistics  of  gun  sys¬ 
tems.  As  such,  investigation  of  the  intricate  behaviour  of  propellant  ignition  and  com¬ 
bustion  on  the  micro-scale  wras  not  a  design  goal.  The  user  is  assumed  to  be  interested 
in  the  macroscopic  behaviour  of  interior  ballistic  arrangements,  with  the  highest  level  of 
accuracy  achievable  with  the  minimum  computational  expense.  A  “fit  for  purpose”  phi¬ 
losophy  dictates  that  complex  physical  models  that  penalise  the  performance  of  the  code 
through  excessive  memory  storage  or  significantly  increased  numerical  operations  should 
only  be  included  when  they  contribute  to  a  significant  increase  in  accuracy  of  the  results 
of  interest. 

It  is  for  this  reason  in  the  context  of  Casbar,  a  decision  has  been  made  not  to  pursue 
the  gas  phase  reaction  model  further  at  this  stage,  due  to  both  the  increased  complexity 
and  the  lack  of  physical  data  for  real  gun  propellants.  The  most  significant  reason  is  the 
lack  of  physicochemical  data  available  on  specific  solid  propellants,  without  which  accurate 
modelling  cannot  be  undertaken.  Some  thermodynamic  data  have  been  published  in  the 
open  literature  ([23]  for  example),  however  the  number  of  propellant  formulations  covered 
is  limited.  It  is  possible  to  establish  a  database  of  the  required  information  for  propellants 
of  interest  to  DSTO,  however  these  tests  would  be  required  to  be  performed  for  each  newT 
formulation  to  render  the  database  accurate  and  up-to-date.  Kinetic  reaction  rate  data 
are  much  more  elusive,  and  difficult  to  determine.  To  develop  reaction  rates  for  the  solid 
propellants  of  interest  for  implementation  into  such  a  model  would  be  a  significant  effort, 
and  as  such  is  not  currently  being  pursued.  Additionally  the  computational  expense  will 
be  greater  than  simpler  models,  reducing  the  performance  of  the  interior  ballistic  code. 
For  these  reasons  the  solid  phase  reaction  model  has  been  chosen  for  implementation  into 
Casbar. 
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3  Numerical  solution  to  the  solid  propellant 

ignition  problem 

This  section  outlines  the  numerical  methods  of  solution  to  the  solid  propellant  ignition 
models  under  investigation  for  use  in  Casbar.  Three  models  of  differing  complexity  are 
investigated  to  determine  the  most  suitable  model.  The  models  can  be  considered  solid- 
phase,  all  chosen  due  to  their  relative  simplicity  in  comparison  to  gas-phase  combustion 
models  and  their  lack  of  requirement  for  reaction  kinetics  in  relation  to  both  heterogeneous 
and  gas  phase  reaction  models.  Diffusion  from  the  oxidizer  into  the  propellant  at  the 
burning  surface  or  within  the  gas  phase  is  not  considered  to  be  a  rate-affecting  process. 
In  addition,  only  one  dimensional  models  are  considered.  Previously  it  has  been  stated 
that  the  burning  rate  of  heterogeneous  propellants  is  strongly  dependent  on  the  particle 
size  of  the  oxidizer  [10],  and  in  such  cases  a  one  dimensional  approximation  would  not  be 
appropriate.  The  majority  of  propellants  considered  in  Casbar  are  homogeneous  single- 
and  double-base  gun  propellants,  wTith  little  variation  with  respect  to  composition  within 
the  grain. 

Two  numerical  approaches  to  solve  the  resulting  one  dimensional  solid  phase  ignition 
model  have  been  investigated.  The  integral  method  is  compared  to  a  solution  b3^  finite 
differences,  with  results  analysed  to  determine  which  solution  would  be  implemented  into 
Casbar. 


3.1  Solution  through  integral  methods 

A  relatively  simple  method  for  solving  the  one  dimensional  heat  conduction  equation  within 
the  solid  phase  involves  using  the  integral  method  [17].  Using  this  method,  the  heat 
penetration  into  the  solid  propellant  grain  is  assumed  to  follow  a  prescribed  profile,  and  is 
dependant  on  the  time  history  of  the  associated  heat  flux,  and  the  instantaneous  heat  flux 
onto  the  propellant  surface.  Once  the  propellant  surface  temperature  reaches  the  defined 
ignition  temperature,  the  grain  is  assumed  to  ignite  and  allowed  to  burn  following  Vieille’s 
pressure-dependant  burning  law.  This  method  is  the  most  commonly  used  in  interior 
ballistic  codes  (eg.  [24],  [25],  [26]  and  [27]). 

The  one  dimensional  heat  conduction  equation  is  significantly  reduced  in  complexity  to 
facilitate  the  solution.  By  the  assumptions  that  all  solid-phase  thermodynamic  properties 
do  not  vary  with  time  or  spatial  position,  that  solid  phase  reactions  are  negligible,  and 
there  is  no  surface  regression  during  the  ignition  transient,  the  heat  conduction  equation 
reduces  to 


dT  ,  d2T 
pC~dt  ~  k7h? 


(12) 


The  method  developed  in  [17]  is  reproduced  here.  By  defining  the  value  <5  (t)  as  the 
time  dependent  penetration  distance  of  the  thermal  wave  front,  the  heat-balance  integral 
can  be  obtained  by  multiplying  the  one-dimensional  heat  conduction  equation  by  dx,  and 
integrating  from  0  to  6, 
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where 


—  (0  +  T0 06)  —  a 


dr  , .  .  dT  .  ; 


(13) 


rm 

6=  /  Tdx  (14) 

Jo 

To  solve  the  resulting  ordinary  differential  equation,  the  temperature  profile  within  the 
solid  propellant  is  assumed  to  take  a  defined  form.  Typically  this  is  of  polynomial  form, 
with  a  cubic  approximation  being  employed  by  [24].  In  this  case  the  temperature  within 
the  solid  profile  at  time  t  is  of  the  form 


T  =  Ax3  +  Bx2  +  Cx  +  D  (15) 

where  A,  B,  C  and  D  are  coefficients  that  may  depend  on  t.  In  order  to  find  expressions 
for  the  coefficients  of  the  polynomial,  boundary  and  initial  conditions  are  required  to  close 
the  system  of  equations. 

As  the  formulation  assumes  a  third-order  profile  to  the  polynomial,  an  additional 
boundary  condition  is  required  to  close  the  system.  This  is  typically  taken  as 

d2T 

Ox*  '*  “  °  (16) 

This  boundary  condition  is  generally  referred  to  as  the  ‘smoothing  function’  as  it  has 
the  effect  of  transitioning  the  temperature  to  the  ambient  conditions  smoothly  [17].  After  a 
suitable  amount  of  algebraic  manipulation  of  the  simultaneous  equations,  the  temperature 
profile  wuthin  the  solid  at  time  t  is  determined  by 


T  (x)  =  J  “  x)3 

where  the  penetration  depth  6  is  given  b3^ 


6  =  \A2ck 


U(*) 


q  (t)  dt 


i/2 


By  setting  x  =  0,  the  surface  temperature  can  be  obtained, 


(17) 


(18) 


Ts  —  Tq  +  \pJjza 


1/2 


q  (■ t )  dt 


(19) 


The  term  q  ( t )  dt  can  be  thought  of  as  the  total  energy  transferred  to  the  propellant 
from  the  moment  of  heat  application  until  time  t  and  is  solved  via  suitable  numerical  tech¬ 
niques,  such  as  the  trapezoidal  rule.  Propellant  ignition  is  assumed  to  have  occurred  when 
the  surface  temperature  reaches  the  (experimentally  determined)  ignition  temperature. 
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3.2  Solution  through  finite  difference  approximation 

Another  approach  to  solving  the  one-dimensional  heat  conduction  equation  for  the  solid 
phase  employs  the  use  of  a  finite  difference  formulation  to  resolve  the  temperature  profile 
within  the  solid  [28].  The  simplified  one  dimensional  heat  conduction  equation  is  again 
used,  with  the  same  bounda^  conditions  applied.  Representing  the  spatial  derivative  with 
a  centred  finite  difference  approximation  at  time  n  and  position  i,  substituting 


d2r  _  +  T^ 

dx2  Ax2 

into  the  solid  phase  conservation  of  energy  equation  gives 


(20) 


dT=  I?+1-21?  +  T?_1 
dt  a  Ax2 


(21) 


The  cell  temperature  at  the  next  time  step  can  then  be  solved  numerically  by  either 
a  first  order  Euler  approximation,  or  a  more  accurate  higher  order  Runge-Kutta  method 
[22],  Using  the  second-order  Runge-Kutta  method,  the  temperature  within  the  solid  at 
the  following  time  step  is  given  by 


w'here 


1  _ 


T T  +  At.f 


(22) 


f  (T?)  =  a 


rpn 

U+ 1 


-  2  Tn 


+  Tn 


i—  1 


Ax2 


(23) 


By  employing  appropriate  boundary  conditions,  the  temperature  profile  within  the 
solid  can  be  established  by  iterating  forward  in  time  from  the  initial  temperature  distribu¬ 
tion,  Too.  Ignition  is  again  established  when  the  surface  temperature  reaches  the  ignition 
temperature. 


3.3  Solution  through  finite  difference  approximation  with  sur¬ 
face  regression 

If  the  one-dimensional  heat  conduction  equation  is  extended  to  include  an  allowance  for 
the  rate  of  energy  convection  from  the  regressing  propellant  surface,  the  conservation  of 
energy  equation  is  then  reformulated  as  the  standard  convection-diffusion  equation, 


dT  dT  ,  d2T 

pcm+pfcdi=k9^ 


(24) 


Again  the  finite  difference  approach  is  employed  to  solve  the  system  of  equations. 
A  suitably  robust  scheme  is  required,  as  throughout  the  ignition  transient  the  problem 
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can  quickly  change  from  diffusion  dominated  to  convection  dominated  causing  significant 
stability  issues  [29].  One  approach  is  to  use  a  forward  in  time  centred  in  space  scheme.  By 
using  a  forward  difference  in  time, 


dT  _  Tj+i  —  T™ 
dt  At 


(25) 


and  centered  in  space, 


d2T  27*.1  -237*  +  77L1 
a— -r  =  a — 


dx 2 


Ax2 


(26) 


dT 

r— —  =  r 
ox 


rpn  _  rpn 

1i+ 1  1i- 1 

2Ax 


(27) 


Combining  into  the  original  PDE,  gives  the  equation  for  the  forward  time  step, 


rpn+l  _  rpn  i  ^  A  t  r  rpn  o  rpn  .  rpn  \  I  ^  ^  ( rTn  rTn  \ 

~  1i  +  Ax2  (J*+l  Z1i  9  A/r  \  *+l 


2Ax 


(28) 


This  scheme  has  been  shown  to  be  unreliable  in  convection  dominated  problems  [29], 
producing  spurious  oscillations  when 


r  Ax 

-  >  2 

a 


However  it  is  locally  stable  when  the  time  step  is  adjusted  such  that  [29] 


(29) 


At  <  min 


Ax2  2a  j 
2a  ’  f2  J 


(30) 


A  suitable  choice  for  time  and  space  discretisation  size  that  avoids  stability  issues  and 
solution  inaccuracies  [29]  is  given  by 


Ax2 

A«<— .  (31) 

where  Ax  also  satisfies 

Ax  <  —  (32) 

~1f 

As  with  the  previous  models,  ignition  is  considered  to  have  occurred  when  the  surface 
temperature  reaches  the  ignition  temperature. 
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4  Solid  phase  ignition  model  solution  technique 

comparison 

In  order  to  develop  a  greater  understanding  of  the  behaviour  of  the  ignition  and  heat 
transfer  models,  a  number  of  synthetic  test  cases  have  been  developed.  Firstly  the  integral 
method  is  compared  to  the  finite  difference  approximation  for  a  number  of  heat  flux  sce¬ 
narios.  A  comparison  is  appropriate,  as  the  same  simplifying  assumptions  are  employed  in 
both.  Both  are  solid  phase  models  that  neglect  the  regression  of  the  solid  surface.  As  such 
it  is  a  comparison  of  the  numerical  solution  technique,  rather  than  the  specific  ignition 
model.  This  therefore  allows  direct  comparison  of  the  numerical  schemes,  possibly  high¬ 
lighting  any  deficiencies  in  the  modelling  behaviour.  The  eight  different  scenarios,  while 
contrived,  will  test  the  models  response  to  heat  flux  input  in  a  variety  of  fashions.  The  first 
four  cases  -  constant  heat  flux,  ramped  heat  flux,  sinusoidal  heat  flux  and  ramped  heat 
flux  -  do  not  model  the  heat  transfer  within  the  packed  propellant  bed.  The  sinusoidal 
heat  flux  and  ramped  heat  flux  scenarios  have  been  included  to  investigate  the  effect  of  a 
varying  flux  on  the  temperature  profiles. 

After  the  simple  solid-phase  reaction  models  are  compared,  the  effect  of  surface  re¬ 
gression  rate,  r,  on  the  ignition  time  of  solid  propellants  is  investigated.  This  is  purely  a 
qualitative  assessment  of  the  effect  of  surface  regression  on  the  ignition  delay  of  a  propel¬ 
lant  grain.  A  range  of  different  values  for  the  regression  rate  are  compared  to  gauge  its 
importance  in  the  transient  ignition  process.  The  propellant  properties  employed  in  this 
study  are  shown  in  Table  1. 


Initial  propellant  temperature,  To 

298  K 

Constant  gas  temperature,  Tg 

2000  K 

Propellant  density,  ps 

1577.8  kg/m3 

Gas  density,  pg 

1.275  kg/m3 

Gas  viscosity,  pg 

1.78  xlO-5  Pa-s 

Gas  ratio  of  specific  heats,  7 

1.4 

Propellant  specific  heat  capacity,  cp 

1550  kJ/kg.I< 

Propellant  thermal  conductivity,  k 

0.31  w/m.K 

Propellant  emissivity,  e 

0.7 

Propellant  ignition  temperature,  Tign 

400  K 

Table  1:  Propellant  properties. 


4.1  Finite  difference  grid  spacing  sensitivity  study 

The  sensitivity  of  the  numerical  solution  to  the  finite  difference  approximation  was  inves¬ 
tigated  to  ensure  solution  independence  from  spatial  discretization.  Three  grid  spacings 
were  investigated,  representing  coarse,  medium  and  fine  grids.  Grid  spacings  of  10  p/rri. 
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lfim  and  0.1  /.im  were  used  to  investigate  the  effect  on  the  surface  temperature  of  the  pro¬ 
pellant  under  the  influence  of  a  constant  heat  flux.  A  plot  of  surface  temperature  versus 
time  is  shown  in  Figure  2,  demonstrating  the  solutions’  dependence  on  grid  spacing.  The 
coarse  grid  shows  a  significant  deviation  from  the  medium  and  fine  grids.  There  is  minimal 
difference  between  the  medium  and  fine  grids.  By  using  the  medium  grid,  a  significant 
saving  in  computing  processing  time  can  be  achieved,  and  hence  has  been  used  as  the  grid 
spacing  for  the  mesh  convergence  quality  assessment. 


Figure  2:  Surface  temperature  sensitivity  to  grid  spacing. 


The  quality  of  spatial  convergence  has  been  quantified  b3^  employing  the  technique 
outlined  in  [30].  By  refining  the  grid  by  a  factor  of  2  around  the  refined  mesh  (Ax  = 
1  x  10_06m),  the  order  of  convergence  can  be  obtained  through 

P  =  ln  (jffZ ~jf)  /In(r) 

where  r  is  the  (constant)  refinement  ratio  and  /  is  the  time  to  ignition  at  the  coarse, 
medium  and  fine  grid  resolutions  (subscripts  3,  2  and  1,  respectively).  By  substituting  in 
the  values  obtained  from  the  numerical  simulations,  the  order  of  convergence  is  calculated 
to  be  1.005.  This  is  somewhat  low  for  a  second  order  in  space  and  time  scheme.  The  grid 
convergence  index  (CGI)  can  be  calculated  as 


CGI  =  Fs 


(A~/i)//i 
(■ rP  -  1) 
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where  Fs  is  a  safety  factor,  taken  as  1.25  when  three  or  more  mesh  spacings  are  em¬ 
ployed.  For  the  medium  and  fine  grids,  the  grid  convergence  index  is  1.53%. 

A  temporal  discretization  study  was  also  performed  for  both  techniques,  and  it  was 
found  that  reducing  the  temporal  discretization  further  from  the  stability  limits  provided 
no  additional  accuracy  to  the  solution. 


4.2  Constant  heat  flux 

The  first  case  investigated  involves  the  comparison  of  the  two  numerical  techniques  under 
the  influence  of  a  constant  heat  flux.  Although  this  is  an  unrealistic  scenario  in  practical 
applications,  it  allows  the  comparison  of  the  two  techniques  with  an  analytical  solution  to 
the  heat  conduction  equation.  The  heat  flux  condition  at  the  boundary  is  given  by 

dT 

k—\x=G  =  q  (33) 

where  q  is  a  constant.  The  propellant  is  assumed  to  have  ignited  when  the  surface 
temperature  reaches  Tign  =400K.  The  analytical  solution  is  given  as  [17] 

Ts  =  T0  +  \Jkj-Kat  •  q  (34) 


Figure  3:  Surface  temperature  time  history,  constant  heat  flux. 

Figure  3  shows  a  plot  of  the  solid  propellant  surface  temperature.  The  plot  shows 
that  the  finite  difference  solution  is  able  to  closely  match  the  analytical  solution,  while 
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the  integral  method  diverges  from  the  actual  solution  as  time  progresses.  The  cause  for 
this  can  be  seen  by  comparing  the  integral  method  cubic  approximation  (Eq.  19)  and  the 
analytical  solution  (Eq.  34).  They  differ  by  approximately  2.3%  at  each  time  step,  which 
increases  the  absolute  error  as  time  progresses. 


4.3  Constant  heat  transfer  coefficient 

The  second  case  investigates  the  individual  models’  response  to  a  temperature  dependent 
heat  flux  at  the  boundary,  such  that 

dT 

k-\x=0  =  h{Tg-Ts(t))  (35) 

where  h  and  Tg  are  held  constant.  The  propellants  are  assumed  to  have  ignited  when 
the  surface  temperature  reaches  400  K.  The  results  of  the  comparison  between  the  cubic 
approximation  and  the  finite  difference  formulation  are  shown  in  Figure  5.  It  can  be  seen 
that  both  curves  have  the  same  general  form,  however  the  finite  difference  approximation 
rises  to  the  ignition  temperature  faster  than  the  cubic  approximation. 


Figure  4 ■  Surface  temperature  time  history,  temperature  dependent  heat  flux. 


4.4  Ramped  heat  heat  transfer  coefficient 

The  third  case  investigates  a  time-dependant  linear  increase  in  heat  transfer  coefficient, 
such  that 
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dT 

k-?^\x=o  =  h{Tg-  Ts(t )) 


(36) 


where 


h  =  f  (t)  =  hconst  ■  t  (37) 

The  results  for  both  the  cubic  approximation  and  the  finite  difference  formulation  are 
shown  in  Figure  5. 


Figure  5:  Surface  temperature  time  history,  ramped  heat  flux. 

The  two  approximations  appear  to  give  similar  results  for  the  ignition  time.  As  in  the 
previous  two  cases,  the  finite  difference  approximation  reaches  the  ignition  temperature  at 
a  faster  rate  than  the  cubic  approximation,  however  the  lag  is  significantly  reduced. 


4.5  Varying  heat  transfer  coefficient 

The  fourth  case  investigated  is  that  of  a  varying  heat  flux  coefficient  with  a  constant 
gas  temperature.  This  case  was  devised  in  order  to  scrutinise  the  behaviour  of  the  cubic 
approximation  against  varying  heat  flux.  The  boundary  condition  is  therefore  given  as 

dT 

k-\x=0  =  h(Tg-Ts(t))  (38) 

where 
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h  =  fit) 


hconst  + 


hconst  sin  (7Tt/180s) 
2 


(39) 


Time,  ms 


Figure  6:  Surface  temperature  time  history,  varying  heat  flux. 


A  plot  of  propellant  surface  temperature  against  time  shows  the  results  from  the  two 
approximations  in  Figure  6.  In  this  case  it  is  apparent  that  the  integral  method  gives 
a  faster  ignition  time  than  the  finite  difference  approximation.  The  integral  method  is 
more  sensitive  to  the  sinusoidal  heat  flux,  which  can  be  seen  in  the  greater  amplitude  of 
the  oscillations  about  the  average  temperature  increase.  The  integral  method  works  by 
assuming  a  temperature  profile  within  the  solid,  based  only  on  the  thermal  properties  of 
the  material,  the  instantaneous  heat  flux  and  the  internal  energy  of  the  solid.  It  does 
not,  therefore,  allow  for  the  lag  of  the  thermal  wave  within  the  solid  due  to  a  decrease  in 
incident  heat  flux.  This  is  an  obvious  shortcoming  of  the  method,  but  the  final  results  are 
not  significantly  affected  in  this  case,  with  the  ignition  time  decreased  by  0.1  ms. 

From  Figure  6  it  is  clear  that  the  difference  in  the  ignition  time  between  the  two 
methods  is  dependant  on  the  choice  of  Tign.  This  could  present  an  issue,  as  at  times 
the  two  methods  may  differ  by  up  to  0.3  ms.  Another  point  that  should  be  noted  is  the 
difference  in  the  peak  value  of  the  oscillations  for  the  solution  as  time  progresses.  This 
implies  that  for  some  ignition  temperatures,  one  method  may  indicate  the  propellant  is 
ignited  while  the  other  method  lags  until  the  next  heat  flux  wave  peak.  This  could  result 
in  significantly  different  ignition  time  estimations  between  the  two  methods. 
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4.6  Step  function  heat  transfer  coefficient 

To  further  test  the  response  of  both  methods  to  highly  variable  heat  fluxes,  a  step  function 
heat  flux  profile  was  employed.  By  inspection  of  the  surface  temperature  profile  equation 
obtained  from  the  integral  method,  the  surface  temperature  will  be  equal  to  the  initial 
temperature  when  the  instantaneous  heat  flux  is  zero.  It  is  therefore  obvious  that  an 
incorrect  surface  temperature  will  result  when  the  heat  flux  is  zero.  It  is  of  interest  however 
to  see  if  the  method  is  able  to  recover  to  the  correct  surface  temperature  when  the  heat  flux 
is  reapplied  based  on  the  total  internal  energy  of  the  propellant.  The  heat  flux  boundary 
condition  is  given  as 


where 


dT 

^lkc^x=0  =  h  ^ 9 


Ts(t)) 


h  =  f  (t)  =  hconstant  '  6  (t)  ; 


15  (t) 


0 

1 

0 


0  <  t  <  0.5 ms 
0.5  <  t  <  1.0ms 
1.0  <  t  <  1.5 ms 
etc 


(40) 


(41) 


Figure  7:  Surface  temperature  time  history,  step  function  heat  flux. 

The  results  from  the  simulations  are  shown  in  Figure  7.  The  finite  difference  approxi¬ 
mation  shows  an  initial  increase  in  surface  temperature,  which  decreases  slowly  when  the 
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heat  flux  is  switched  off.  As  the  finite  difference  technique  solves  the  temperature  profile 
within  the  solid,  heat  diffusion  away  from  the  surface  can  be  calculated.  The  behaviour  of 
the  integral  method  to  this  (quite  contrived)  heat  flux  is  significantly  different.  The  two 
methods  follow'  a  very  similar  profile  upon  the  initial  application  of  heat,  but  when  the 
flux  is  switched  off  the  surface  temperature  approximated  by  the  integral  method  drops 
immediately  to  zero.  Upon  the  re-application  of  heat,  the  surface  temperature  given  by 
the  cubic  approximation  jumps  up  to  its  previous  position.  The  temperature  profile  within 
the  solid  is  derived  from  the  stored  internal  energy  within  the  solid,  and  the  instantaneous 
heat  flux  to  the  surface.  The  integral  method  over-estimates  surface  temperature  at  this 
point  (as  compared  to  the  finite  difference  equation),  which  is  attributed  to  the  assumption 
of  a  cubic  profile  no  longer  being  valid.  While  there  is  no  heat  flux  to  the  solid  propellant 
surface,  heat  is  still  conducted  away  from  the  hot  solid  regions  to  cooler  regions  further 
from  the  propellant  surface.  The  temperature  profile  is  therefore  no  longer  cubic,  and  a  key 
assumption  for  the  formulation  of  the  cubic  approximation  is  no  longer  valid.  The  cubic 
approximation  gets  further  away  from  the  finite  difference  approximation  as  the  number  of 
flux  intervals  increases.  For  the  specified  Tign,  by  chance,  the  final  ignition  time  predicted 
by  the  two  methods  is  very  similar. 


4.7  Packed  bed  heat  transfer  coefficient,  varying  velocity  and 
constant  temperature 

The  next  step  in  the  analysis  involved  implementation  of  a  packed  bed  heat  transfer  approx¬ 
imation.  This  wras  done  in  order  to  assess  both  models’  behaviour  under  the  application  of 
more  complex  and  realistic  heat  flux  scenarios.  The  heat  transfer  is  a  function  of  a  number 
of  physical  processes  all  acting  to  transfer  heat  energy  to  the  solid  propellant.  Due  to  the 
complexity  of  the  heat  transfer  mechanisms  present,  no  one  complete  method  exists  that 
accurately  accounts  for  all  processes.  Rather,  empirical  formulations  are  typically  used  to 
approximate  the  amount  of  heat  transferred  to  (and  from)  the  solid  propellant  and  the 
surroundings.  This  section  outlines  some  of  the  models  that  have  been  previously  used  in 
interior  ballistic  and  propellant  combustion  modelling,  with  the  aims  of  selecting  the  most 
suitable  model  for  use  in  Casbar. 

Gough  [31]  employs  an  empirical  approximation  to  the  heat  transfer,  assuming  radiative 
and  convective  mechanisms,  based  on  fluidized  bed  theory  [32],  The  heat  transferred  to 
the  solid  propellant  grain  is  a  function  of  temperature  and  is  given  by 

q  =  ( hconv  +  hrad )  (Tg  -  Ts)  (42) 

The  convective  and  radiative  film  coefficients  are  determined  by 

Kad  =  ea  (Tg  +  Ts)  (T2g  +  T2)  (43) 


and 


h 


conv 


Nu  •  kf 
Dp 


(44) 
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In  the  above  expressions  e  is  the  emissivity  of  the  propellant /gas  interface,  a  is  the 
Stefan  Boltzmann  constant,  Nu  is  the  Nusselt  number,  kf  is  the  thermal  conductivity  of 
the  gas  evaluated  at  the  film  temperature  and  Dp  is  the  effective  grain  diameter,  which  is 
given  by 


DP  = 


6  Vp 
SP 


The  Nusselt  number  is  determined  by 


Nu  =  0APr1/3Re2/3 


(45) 


(46) 


where 


Rcs  —  pf\ug  us\Dp/ pj 


(47) 


and 


Pr 


Cs^f 

kf 


(48) 


In  the  above  expression  n f  and  kf  are  the  viscosity  and  thermal  conductivity  of  the  gas 
evaluated  at  the  film  temperature  and  cs  is  the  specific  heat  capacity  of  the  solid  phase. 
Lowe  [22]  also  employs  an  empirical  approximation  for  the  heat  transfer,  however  the  heat 
transfer  is  assumed  to  occur  only  through  radiation  and  conduction,  with  convective  effects 
neglected.  This  inaccuracy  is  acknowledged,  with  experiments  highlighting  the  need  for  a 
scaling  factor  to  be  introduced  to  obtained  realistic  solid  temperature  profiles.  Identical  to 
[31],  the  heat  flux  due  to  radiation  is  given  by 


Qrad  =  ve  (T4  -  T4)  (49) 

The  conductive  heat  transfer  differs  however,  and  is  given  by  the  empirical  relation 

qcond  =  kfcsT°-25  7^(Tg-Ts)  (50) 

up 

No  details  as  to  the  origin  or  background  theory  of  this  expression  were  provided  in 

[22]- 

Both  of  the  above  approaches  to  modelling  the  heat  transfer  within  the  packed  propel¬ 
lant  bed  neglect  all  mechanisms  besides  convective  and  radiative  heating.  Unfortunately 
this  approach  neglects  the  heat  transfer  due  to  non-gaseous  igniter  products.  These  prod¬ 
ucts  can  form  a  significant  portion  of  the  igniter  gases,  and  as  such  the  heat  transfer  can 
be  under  predicted.  Particle-particle  interaction  models  can  be  employed  to  model  this 
aspect  of  heat  transfer  [33],  however  this  approach  is  not  currently  being  pursued  within 
the  Casbar  framework. 
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The  first  numerical  simulation  involved  examining  the  behaviour  under  the  influence  of 
a  varying  gas  velocity,  with  a  constant  gas  temperature.  The  heat  flux  boundary  condition 
is  therefore 

dT 

k-\x=0  =  h(Tg-Ts(t))  (51) 

where 


h 


hconv  T  hrad 


and 


Nuk 


hrad  —  ec  (' Tg  +  Ts)  (Tg  +  T2) 


(52) 


(53) 

(54) 


Figure  8:  Surface  temperature  time  history,  packed  bed  heat  transfer  approximation  with 
varying  velocity. 


The  results  shown  in  Figure  8  indicate  that  by  varying  the  gas  velocity,  little  influence  is 
made  on  the  surface  temperature  of  the  propellants  calculated  by  the  two  methods.  This  is 
simply  because  convective  heat  transfer  plays  a  much  less  significant  role  in  the  propellant 
ignition  case  being  analysed. 
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4.8  Packed  bed  heat  transfer  coefficient,  constant  velocity 
and  varying  temperature 

The  final  case  being  examined  is  again  using  the  packed  bed  heat  transfer  approximations, 
however  this  time  with  a  varying  gas  temperature.  In  this  test  case,  the  gas  temperature  is 
varied  by  a  sinusoidal  function,  giving  oscillations  about  a  defined  mean  temperature.  As 
radiation  plays  a  much  more  significant  role  in  the  heat  transfer  to  the  solid  propellant  at 
high  temperatures  due  to  the  quartic  dependencj^,  both  approximations  will  be  significantly 
tested  by  this  condition. 


Figure  9:  Surface  temperature  time  history,  packed  bed  heat  transfer  approximation  with 
varying  gas  temperature. 


Figure  9  shows  the  response  of  the  two  models  to  the  applied  heat  flux.  It  can  be  clearly 
seen  that  the  integral  is  significantly  affected  by  the  varying  heat  flux.  This  is  again  due 
to  the  sensitivity  of  the  method  to  the  instantaneous  heat  flux.  The  final  ignition  times 
determined  by  both  methods  does  not  vary  significantly  between  the  two  methods. 

4.9  Effect  of  surface  regression  on  ignition  time 

To  understand  the  effect  of  surface  regression  on  the  ignition  time  of  the  propellant,  the 
finite  difference  case  discussed  in  Section  3.3  with  surface  regression  has  been  investigated. 
Figure  10  shows  surface  temperature  history  plots  for  various  values  of  f,  as  well  as  the 
original  non-regressing  case.  As  the  regression  rate  increases,  the  time  to  ignition  of  the 
solid  propellant  also  increases.  This  can  be  attributed  to  the  convection  term  within  the 
energy  equation. 
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Figure  10:  Effect  of  varying  regression  rate  on  propellant  ignition  time. 


4.10  Numerical  scheme  performance  comparison 

A  comparison  of  the  integral  and  finite  difference  methods  for  the  propellant  ignition  prob¬ 
lem  was  performed  on  the  constant  heat  transfer  scenario.  As  this  is  the  simplest  heat 
transfer  problem,  it  will  provide  somewhat  of  a  ‘best  case  scenario’  in  terms  of  computa¬ 
tional  efficiency.  This  calculation,  which  can  be  considered  typical  for  a  solid  propellant 
ignition  calculation  within  Casbar,  would  highlight  the  efficiency  differences  of  the  two 
methods.  Both  numerical  schemes  were  employed  with  the  same  initial  and  boundary  con¬ 
ditions,  with  the  processing  time  recorded.  For  the  integral  method,  a  single  calculation 
of  the  surface  temperature  would  take  0.574  ms,  as  opposed  to  2.69  ms  for  the  finite  dif¬ 
ference  approach.  This  represents  a  four-fold  increase  in  the  computational  expense  of  the 
finite  difference  method  over  the  integral  method  based  on  calculation  time  alone.  Consid¬ 
ering  that  this  calculation  is  required  to  be  performed  in  each  individual  cell  at  every  time 
step,  a  significant  computational  penalty  is  associated  with  employing  the  finite  difference 
approach. 


5  Summary  and  conclusion 

A  wide  range  of  approaches  to  solid  propellant  ignition  modelling  are  available.  The  physi- 
cal  process  of  ignition  and  combustion  in  solid  propellants  is  complex  and  case- dependant, 
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and  as  such  no  single  model  can  be  considered  appropriate  for  every  problem.  Simplifica¬ 
tions  are  required,  which  may  or  may  not  be  applicable  for  different  scenarios.  The  three 
main  areas  of  ignition  models  (solid-phase,  heterogeneous  and  gas-phase  reactions  models) 
encompass  a  broad  range  of  complexity  and  numerical  efficiency.  It  is  therefore  desirable  to 
choose  a  model  that  is  accurate,  while  being  flexible  for  a  wide  range  of  possible  scenarios 
for  implementation  into  Casbar. 

In  Casbar  it  is  required  that  the  interior  ballistic  propellant  ignition  model  is  adequately 
robust,  while  reducing  the  impact  on  the  simulation  time.  In  all  of  the  investigated  heat 
flux  scenarios,  the  integral  method  was  able  to  adequately  approximate  the  final  ignition 
time  to  within  an  acceptable  level  of  accuracy  (in  comparison  to  the  finite  difference  ap¬ 
proximation).  Situations  involving  highly  variable  heat  fluxes  were  employed  to  test  the 
applicability  of  the  integral  method.  These  scenarios  were  constructed  to  accentuate  the 
variability  of  the  heat  flux,  and  situations  like  this  are  not  expected  in  typical  interior  bal¬ 
listic  simulations.  The  integral  method  is  therefore  considered  an  appropriate  candidate 
for  implementation  in  Casbar. 
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